PL-TR-96-2299 

Environmental  Research  Papers,  No.  1201 


ENHANCEMENTS  TO  THE  CROSS  CORRELATION 
TECHNIQUE  FOR  EXTRAPOLATING 
GEOSTATIONARY  SATELLITE  IMAGERY 


Kenneth  F.  Heideman 
R.  Radburn  Robb,  1  Lt,  USAF 
Donald  C.  Norquist 


27  November  1996 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED. 


19971121  016 

[PTIC  8 


PHILLIPS  LABORATORY 

Directorate  of  Geophysics 

AIR  FORCE  MATERIEL  COMMAND 

HANSCOM  AIR  FORCE  BASE,  MA  01731-3010 


"This  technical  report  has  been  reviewed  and  is  approved  for  publication." 


/(  (ytuU  C- 

DONALD  C.  NORQUISJ  DONALD  A.  CHT^HOLM 

Acting  Chief,  Satellite  Analysis  and  Director,  Atmospheric  Sciences  Division 
Weather  Prediction  Branch 
Atmospheric  Sciences  Division 


This  report  has  been  reviewed  by  the  ESC  Public  Affairs  Office  (PA)  and  is  releasable 
to  the  National  Technical  Information  Service  (NTIS). 


Qualified  requestors  may  obtain  additional  copies  from  the  Defense  Technical 
Information  Center.  All  others  should  apply  to  the  National  Technical  Information 
Service. 


If  your  address  has  changed,  or  if  you  wish  to  be  removed  from  the  mailing  list,  or  if  the 
addressee  is  no  longer  employed  by  your  organization,  please  notify  PL/IM,  29 
Randolph  Street,  Hanscom  AFB,  MA  01731-3010.  This  will  assist  us  in  maintaining  a 
current  maiiing  list. 


Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices  on  a  specific 
document  requires  that  it  be  returned. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response.  Including  the  time  for  reviewing  instructions,  searchir>g  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1216  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  10704-0188),  WasNngton,  DC  20B03. 


1.  AGENCY  USE  ONLY  (Leave  blanks  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

27  November  1996  Final  Report  October  94  to  October  96 


4,  TITLE  AND  SUBTITLE 

Enhancements  to  the  Cross  Correlation  Technique  for  Extrapolating  Geostationary 
Satellite  Imagery 


6.  AUTHOR(S) 

Heideman,  Kenneth  F. 

Robb,  R.  Radbum,  ILt,  USAF 

Norquist,  Donald  C.  _ 


7.  PERFORMING  ORGANIZATION  NAMEIS)  AND  ADDRESS(ES1 

Phillips  Laboratory  (GPAB) 

29  Randolph  Road 

Hanscom  AFB,  MA  01731-3010 


5.  FUNDING  NUMBERS 

PE  62101 
PR  6670 
TA  GT 
WU35 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

PL-TR-96-2299 
ERP,  No.  1201 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS{ES) 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


12a.  DISTRIBUTION  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words} 

This  report  documents  enhancements  made  to  the  cross  correlation  technique  for  extrapolating  geostationary  satellite  imagery 
up  to  five  hours  into  the  future.  Because  it  requires  only  two  geostationary  images  to  utilize  the  technique,  it  satisfies  an 
important  Air  Force  requirement  to  be  able  to  forecast  the  distribution  of  cloud  in  data  restricted  environments.  The 
enhancements  referred  to  above  include  techniques  to  prevent  "terrain  advection,"  reduce  border  contamination,  and  to 
minimize  excessive  smoothing  of  cloud  features.  In  addition,  an  interactive  Graphic  User  Interface  has  been  developed  to 
facilitate  use  of  the  technique  by  operational  weather  forecasters. 


14.  SUBJECT  TERMS 

Cloud 

Satellite  imagery 

Extrapolation 

Cross  correlation 

17.  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 


15.  NUMBER  OF  PAGES 
26 


16.  PRICE  CODE 


19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACTi 
OF  ABSTRACT 


Unclassified 


Standard  Form  298  (Rev.  2-89)  (EG) 

Prescribed  by  ANSI  Std.  239.18 

Designed  using  Perform  Pro,  WHS/DIOR,  Oct  94 


Contents 


1.  INTRODUCTION  1 

2.  ENHANCEMENTS  TO  THE  CROSS  CORRELATION  TECHNIQUE  3 

2.1  Terrain  Advection  3 

2.2  Border  Contamination  8 

2.3  Excessive  Smoothing  of  Cloud  Features  9 

2.4  Development  and  Dissipation  of  Cloud  10 

2.4.1  Climatological  Blending  H 

2.4.2  Harmonic  Analysis  12 

3.  DEVELOPMENT  OF  A  GRAPHIC  USER  INTERFACE  FOR  THE 

CROSS  CORRELATION  TECHNIQUE  13 

3.1  Porting  and  Code  Optimization  14 

3.2  Development  of  a  User-Friendly  Interface  14 

3.3  Development  of  Documentation,  Help,  and  Tutorial  15 

4.  CONCLUDING  REMARKS  15 

REFERENCES  17 

APPENDIX:  HARMONIC  ANALYSIS  METHOD  19 


Illustrations 


1.  Cross  Correlation  Technique  Regions  of  Interest:  Eastern  Asia, 

Eastern!  Mediterranean,  and  Central  and  Northern  South  America  5 

2.  Example  of  Hourly  Clearscene  Image  for  Eastern  Mediterranean  ROI  6 

3.  Example  of  Hourly  Clearscene  Image  for  Eastern  Asia  ROI.  Note  cloud 
in  upper  portion  of  imager,  indicating  that,  for  this  particular  hour, 

a  number  of  pixels  were  never  actually  clear  throughout  the  available 

10-day  climatological  sample  of  geostationary  images  7 


IV 


Acknowledgments 


The  authors  thank  Donald  A.  Chisholm,  Acting  Director  of  the 
Atmospheric  Sciences  Division  (PL/GP),  for  his  guidance  and  suggestions  in 
further  developing  and  testing  the  application  of  the  cross-correlation  technique  to 
forecast  the  time-evolution  of  geostationary  imagery. 


VI 


Enhancements  to  the  Cross  Correlation  Technique 
for  Extrapolating  Geostationary  Satellite  Imagery 

1.  INTRODUCTION 


Accurate  nowcasting  (0-2  hours)  and  very  short-range  forecasting  (up  to  6 
hours)  of  cloud  features  has  long  been  a  top  priority  of  the  Air  Force.  Knowledge 
of  the  distribution  of  cloud  is  essential  for  communication,  training  exercises, 
tactical  deployment  of  military  aircraft,  and  surveillance.  Such  forecasts  must 
often  be  made  with  a  minimum  of  observational  data  to  accommodate  tactical 
situations  in  which  observations  are  largely  unavailable  or  detection  must  be 
prevented. 

Since  1989,  developing,  evaluating,  and  refining  methods  for  data-sparse 
forecasting  has  been  an  important  component  of  the  research  conducted  within  the 
Atmospheric  Sciences  Division  of  Phillips  Laboratory  (PL/GPA).  Early  efforts 
focused  on  identification,  isolation,  and  extrapolation  of  individual  cloud  features 
based  on  infrared  brightness  temperature  contours.  Heideman  et  al  (1990), 
Ruggiero  and  Heideman  (1992),  and  Nehrkom  et  al  (1993)  have  documented  the 
details  and  evaluation  of  several  contour  extrapolation  techniques. 

In  general,  published  results  suggest  that,  with  respect  to  simple  persistence, 
some  success  in  forecasting  the  shape,  movement,  and  evolution  of  cloud  features 
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was  achieved  using  contour  extrapolation,  but  the  disadvantages  were  significant. 
For  example,  none  of  the  techniques  studied  were  designed  to  account  for  the 
natural  tendency  of  individual  cloud  features  to  split  and  merge  with  other  clouds 
over  time.  Moreover,  Heideman  et  al  (1990)  detail  the  cumbersome,  largely 
manual,  user-intensive  contour  selection  process  required  to  process  an  image  and 
produce  a  forecast.  The  computational  and  manual  expense  precluded  researchers 
from  evaluating  a  sufficient  number  of  contour  extrapolation  cases  to  make  a 
definitive  recommendation  to  either  expand  research  on  the  use  of  one  or  more  of 
the  candidate  techniques,  or  to  move  in  a  different  research  direction. 

As  a  result,  in  1992  PL  initiated  a  project  to  test  the  established  candidate 
contour  extrapolation  techniques,  as  well  as  to  investigate  other  possible 
techniques,  pursuant  to  declaring  one  technique  most  suitable  for  further 
development  and  potential  operational  implementation.  Detailed  findings  are 
documented  by  Nehrkom  et  al  (1993).  A  whole-image  extrapolation  technique 
was  developed  in  the  comparisons  study.  This  method  was  based  on  the  statistical 
technique  of  cross  correlation.  It  was  demonstrably  superior  to  the  candidate 
contour  extrapolation  techniques  as  well  as  to  persistence.  The  technique  required 
only  two  consecutive  geostationary  visible  or  infrared  satellite  images,  30  or  60 
minutes  apart,  to  produce  forecast  images  up  to  5  hours  into  the  future.  This 
means  that  the  technique  can  be  used  in  data-restricted  circumstances. 

Infrared  imagery  has  been  used  exclusively  in  this  study  to  take  advantage  of 
its  24-hour  utility.  For  display  purposes,  pixel  infi'ared  brightness  temperatures  are 
routinely  converted  to  grayshade  units  using  a  linear  scale.  Forecasts  were  the 
product  of  displacement  vectors  based  on  differences  in  pixel  position  and 
brightness  between  the  two  observed  images  needed  to  initiate  the  technique. 

Beyond  better  verification  scores,  cross  correlation  had  a  two-fold  advantage 
over  the  other  techniques  considered.  First,  it  implicitly  accounted  for  splitting  and 
merging  cloud  features.  Second,  forecast  products  could  be  looped  like  a  series  of 
actual  geostationary  images.  For  specific  details  of  how  the  cross  correlation 
technique  works,  and  how  geostationary  imagery  must  be  preprocessed  prior  to 
running  the  program,  see  Hamill  and  Nehrkom  (1992)  and  Nehrkom  et  al  (1993). 

Cross  correlation  was  selected  as  the  cloud  extrapolation  technique  vdth  the 
greatest  potential  for  further  development  and  eventual  operational 
implementation.  However,  there  were  obvious  limitations  of  and  problems  with 
the  cross  correlation  technique.  For  example,  in  largely  clear  scenes,  terrain 
features  such  as  lakes  and  rivers  were  often  interpreted  as  cloud  because  their 
infi'ared  brightness  temperatures  were  significantly  different  from  their 
surroundings.  As  a  result,  they  would  often  be  "advected"  downstream  in  forecast 
satellite  images.  This  came  to  be  known  as  the  "terrain  advection”  problem.  In 
addition,  contamination  at  the  edges  of  the  forecast  images  often  became  a 
problem.  Finally,  the  technique  was  unable  to  represent  or  infer  cloud 
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development  or  dissipation.  These  are  among  the  major  concerns  that  have  been 
addressed  in  this  report.  We  will  also  discuss  an  effort  to  develop  an  operational 
user  interface  for  the  cross  correlation  technique. 

2.  ENHANCEMENTS  TO  THE  CROSS  CORRELATION  TECHNIQUE 

Nehrkom  et  al  (1993)  recommended  that  fixture  work  on  the  cross  correlation 
technique  be  focused  on  addressing  several  key  issues:  terrain  advection,  border 
contamination,  the  loss  of  image  contrast  with  time,  and  the  lack  of  cloud 
development  and  dissipation.  The  following  sections  describe  recent  elTorts  to 
solve  these  problems. 

2.1  Terrain  Advection 

As  Hamill  and  Nehrkom  (1992)  explain,  the  cross  correlation  technique 
employs  a  backwards  trajectory  scheme,  whereby  pixel  intensities  upstream  (that 
is,  at  the  origin  of  displacement  vectors)  become  the  forecast  pixel  intensity  values. 
Problems  arise  when  the  upstream  pixels  are  over  bodies  of  water,  which  often 
have  brightness  signatures  substantially  different  than  the  surrounding  terrain.  In 
these  situations,  oceans,  lakes,  and  rivers  are  often  interpreted  as  cloud,  and  the 
result  is  that  they  are  advected  downstream  in  the  forecast  imagery.  Clearly,  it  is 
untenable  to  forecast  that  Lake  Michigan,  for  example,  will  move  to  the  east  within 
a  five-hour  period. 

The  solution,  as  proposed  by  Hamill  and  Nehrkom,  is  to  use  a  nephanalysis 
valid  at  the  time  of  the  most  recent  observed  satellite  image.  This  nephanalysis 
provides  information  on  which  pixels  were  cloudy  and  which  were  clear.  Because 
that  image  is  the  one  fi’om  which  all  forecast  images  are  generated,  a  check  can  be 
made  for  each  pixel  to  determine  whether  its  upstream  counterpart  was  deemed 
clear  or  cloudy  in  the  nephanalysis.  To  avoid  terrain  advection,  upstream  clear 
pixel  intensities  are  simply  not  advected  to  the  forecast  pixel  location;  rather,  a 
clearscene  climatology  pixel  intensity  value  becomes  the  forecast  value.  This  works 
well  in  theory,  but  it  depends,  of  course,  on  the  availability  of  a  nephanalysis  and 
representative  clearscene  mask  for  the  study  area.  Moreover,  the  quality  of  the 
forecast  images  is  obviously  dependent  on  the  accuracy  of  the  nephanalysis.  The 
availability  of  accurate  nephanalysed  data  has  benefitted  from  the  Support  of 
Environmental  Requirements  for  Cloud  Analysis  and  Archive  (SERCAA)  cloud 
analysis  products  (Isaacs,  1993  and  Neu  et  al.,  1994).  Most  of  the  geostationary 
imagery  used  in  this  project  was  made  available  through  the  SERCAA  database  at 
PL.  One  consequence  of  this  has  been  the  60-minute  interval  between  images 
archived  for  SERCAA,  as  compared  to  the  half-hourly  images  used  by  Hamill  and 
Nehrkom  (1992).  As  expected,  the  forecasts  produced  with  the  hourly  images 
generally  have  larger  RMSEs  than  the  half-hourly  images,  but  the  difference  is 
small  thi^ough  2  hours  (roughly  3%,  on  average).  No  direct  comparative  data  exist 
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beyond  2  hours,  as  forecasts  made  using  the  half-hourly  images  were  limited  to  2.5 
hours,  as  opposed  to  5  hours  using  hourly  images. 

The  primary  task  in  Phase  I  of  SERCAA,  already  completed,  was  to  integrate 
multispectral  imagery  from  geostationary  and  polar  orbiting  satellite  platforms. 
Algorithms  were  applied  to  images  from  each  individual  platform  and  then  to  the 
integrated  images  to  determine  the  distribution  of  cloud,  both  in  a  binary  (0  =  no 
cloud,  1  =  cloud  present)  sense  (on  a  pixel  by  pixel  basis),  and  in  terms  of  cloud 
fraction  on  a  l/16th-mesh  basis  (251cm  resolution). 

The  binary  pixel  by  pixel  analyses  have  been  used  in  this  project  to  prevent 
terrain  advection.  SERCAA  nephanalyses  were  generated  hourly  over  10-day 
periods  throughout  1992  and  1993  in  several  regions  of  interest  (ROIs).  The 
Eastern  Mediterranean  and  Eastern  Asian  ROIs  were  used  in  evaluating  and 
refining  the  cross  correlation  technique,  and  the  Central  American  ROI  was  used  in 
subsequent  testing  and  evaluation  (see  Figure  1). 

The  clearscene  climatology  consists  of  a  single  image  for  each  hour  (Figure  2). 
It  was  created  by  identifying  the  minimum  pixel  intensity  values  on  an  hourly  basis 
over  the  10-day  period  composing  each  data  set,  and  assigning  that  grayshade  to 
that  particular  location  in  the  clearscene.  This  is  not  an  ideal  solution,  as  by 
definition  the  clearscene  pixels  represent  the  low-end  extreme  grayshade  values  for 
a  given  hour  and  location  over  a  10-day  period,  and  these  pixels  will  often  appear 
artificially  darker  than  the  corresponding  verification  satellite  imagery.  In  tropical 
areas  (for  example,  the  eastern  Asia  ROI),  cloud  is  often  ubiquitous,  and  the 
minimum  grayshade  for  a  given  hour  over  the  10-day  period  ac^ally  represents 
cloud  for  a  number  of  pixels.  This  results  in  a  slightly  "milky"  appearance  of  the 
clearscene  image  (see  Figure  3  ).  A  possible  solution  for  this  is  to  set  a  maximum 
threshold  for  grayshade  on  the  clearscene  image  that  is  less  than  the  value  likely  to 
be  associated  with  cloud  in  the  area  of  interest. 

Alternatively,  it  is  possible  to  substitute  color-coded  pixels  based  on  geography 
(land,  coastline,  water)  for  clearscene  pbcel  values,  using  a  1/1 6th  mesh  resolution 
hemispheric  database  developed  at  PL  for  use  in  SERCAA.  This  "geography 
option"  is  currently  being  coded  and  will  be  explored  as  a  viable  operational 
alternative  to  the  clearscene  mask.  The  clearscene  problem  notwithstanding,  the 
basic  approach  of  not  advecting  upstream  clear  pixels  is  an  unqualified  success  in 
preventing  terrain  advection.  Testing  of  the  technique  over  the  eastern 
Mediterranean  ROI,  which  contains  portions  of  the  Mediterranean  and  Red  Seas, 
and  Lake  Nasser  (easily  distinguishable  large  bodies  of  water  prone  to  spurious 
advection  if  a  nephanalysis  is  not  used),  confirms  that  the  approach  of  preventing 
advection  of  upstream  clear  pixels  effectively  eliminates  the  problem.  The  added 
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Figure  1;  Cross  Correlation  technique  Regions 
of  Interest:  Eastern  Asia,  Eastern  Mediterranean, 
and  Central  and  Northern  South  America. 


Figure  2:  Example  of  hourly  clearscene  image  for  Eastern  Mediterranean  ROI. 


Figure  3:  Example  of  hourly  clearscene  image  for  Eastern  Asia  ROI.  Note  cloud 
in  upper  portion  of  imager,  indicating  that,  for  this  particular  hour,  a  number  of 
pixels  were  never  actually  clear  throughout  the  available  ten-day  climatological 
sample  of  geostationary  images. 


computational  cost  is  negligible  when  compared  to  the  benefit  of  ensuring  an 
anchored  background  under  moving  clouds  in  forecast  animation  loops. 

2.2  Border  Contamination 

Inaccuracies  and  distortions  along  the  boundary  of  any  limited-area  domain  are 
to  be  expected.  Advection  of  pixel  gray  shade  values  based  on  a  recent  history  of 
their  positions  is  the  essence  of  the  cross  correlation  technique.  There  is  no  such 
information  outside  the  image  domain,  and  advection  right  at  the  boundary  is  zero. 
Often,  the  result  is  distorted  cloud  features  at  the  edge  of  the  domain  that  appear 
blurred  and/or  fail  to  move  in  synch  with  the  rest  of  the  forecast  images  when 
looped.  Cloud  features  that  existed  just  outside  the  image  domain  at  the  time  of 
the  latest  observed  image  cannot  be  incorporated  in  forecast  images. 

Our  solution  to  the  border  contamination  problem  was  to  enlarge  the  domain 
over  which  the  technique  was  applied  but  to  display  only  the  inner  domain  of  the 
original  size.  Prior  to  addressing  this  problem,  we  had  uniformly  worked  with  256 
X  256  pixel  infrared  images,  with  resolutions  varying  between  4  km  for 
METEOSAT  (Europe  and  environs)  and  GMS  (greater  Asia),  and  8-km  for 
GrOES-7  (Americas).  The  question  of  how  large  the  additional  buffer  should  be 
was  important,  as  the  computational  expense  of  increasing  the  image  domain  was 
non-trivial.  The  primary  consideration  was  that  it  be  large  enough  that 
contamination  at  the  edge  of  the  buffer  domain  in  the  first  forecast  image  not 
advect  into  the  inner  domain  during  the  5-hour  forecast  period.  Assuming  an 
average  advective  speed  of  80  km/hr  and  basing  our  calculations  on  the  best 
available  resolution  of  4  km,  we  calculated  a  displacement  of  100  pixels  (400  km) 
in  5  hours.  That  was  the  minimum  desired  width  of  the  buffer  zone.  We 
conservatively  settled  on  a  buffer  zone  128  pixels  wide  to  add  to  the  original  256  x 
256.  As  a  result,  the  working  domain  became  384  x  384  pixels,  thus  providing  a 
buffer  of  64  pixels  on  each  side  of  the  displayed  256  x  256  image. 

Testing  suggested  that  the  approach  is  very  effective,  with  whole-image  root 
mean  square  errors  (RMSE)  reduced  an  average  of  15-20  percent;  the  subjective 
aesthetic  improvement  seemed  far  greater.  The  best  improvement  was  for  those 
images  with  a  great  deal  of  cloud  along  the  edge  or  just  outside  of  the  observed 
256  X  256  images. 

The  cost  of  the  improvement,  however,  has  been  a  20-25  percent  increase  in 
computer  time,  bringing  the  run  time  to  nearly  25  minutes  on  a  VAX4000 
workstation  for  a  5-hour  forecast.  Reducing  the  width  of  the  buffer  would 
certainly  help  in  this  regard,  but  the  resulting  increased  probability  that  distorted 
cloud  would  find  its  way  into  the  later  forecast  images  makes  that  a  questionable 
strategy.  A  better  solution  would  be  to  optimize  the  code  to  speed  up  execution  of 
the  program  (see  Section  3  a)  and  to  use  a  more  powerful  computer  to  execute  the 
cross  correlation  technique. 
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2.3  Excessive  Smoothing  of  Cloud  Features 

Of  the  four  basic  problems  with  the  cross  correlation  technique  mentioned  at 
the  beginning  of  this  section,  we  feel  that  the  slight  loss  of  contrast  over  time  in  the 
forecast  images  is  the  least  serious.  However,  it  was  certainly  worth  solving  in  our 
attempt  to  produce  the  most  realistic  looking  forecast  images  possible. 

The  smoothing  problem  results  from  excessive  smoothing  and  averaging  in  the 
extrapolation  process.  The  key  to  remedying  it  lies  in  determining  at  what  points 
in  the  process  of  generating  the  forecast  images,  if  any,  averaging  is  unnecessary 
and  can  be  eliminated  without  adversely  affecting  the  end  product.  Averaging  is 
required  in  calculating  the  cross  correlation  displacement  vector  field,  running  the 
quality  control  filter,  and  performing  the  objective  analysis  to  assign  displacement 
vectors  at  every  pixel  (Hamill  and  Nehrkom,  1992).  However,  there  was  some 
flexibility  in  two  other  components  of  the  program  that  applied  varying  degrees  of 
smoothing  and  averaging  to  the  forecast  satellite  images:  the  use  of  bilinear 
interpolation  to  produce  forecast  grayshade  values  and  the  number  of  times  it  is 
used  when  tracing  a  grayshade  back  to  its  original  location  (a  necessary  step  in 
determining  the  value  of  the  pixel  in  the  next  forecast  frame). 

When  using  the  backwards  trajectory  scheme  described  in  Hamill  and 
Nehrkom  (1992)  to  obtain  forecast  grayshade  values,  displacement  vectors  are 
followed  back  to  their  origin  in  an  earlier  frame,  which  may  or  may  not  lie  directly 
on  a  pixel.  Bilinear  interpolation,  which  is  a  weighted  average  of  the  surrounding 
four  pixel  values,  has  been  used  to  compute  the  forecast  pixel  value  in  such  cases. 
This  is  certainly  a  reasonable  approach,  but  may  result  in  a  subtle  smoothing  of 
cloud  features  on  forecast  imagery  that  increases  with  increasing  forecast  interval. 
As  Hamill  and  Nehrkom  point  out,  one  alternative  would  be  to  round  the 
trajectory  calculation  to  the  nearest  pixel  instead  of  using  bilinear  interpolation; 
this  would  result  in  sharper  images,  but  likely  at  the  expense  of  accuracy. 

To  test  this  notion,  eight  cross  correlation  cases  were  mn  with  and  without 
bilinear  interpolation,  using  30  and  60  minute  intervals  between  observed  images. 
RMSEs  were  computed  separately  for  each  of  the  five  forecast  images  and  each  of 
these  were  averaged  over  the  eight  cases.  Results  are  shown  in  Table  1.  It  is 
apparent  that  use  of  bilinear  interpolation  reduced  RMSE  regardless  of  the  time 
interval  between  images  and  at  all  forecast  times,  on  average,  by  as  much  as  12.3 
percent  for  a  2.5-hour  forecast  for  the  30-minute  update  images,  and  by  as  much 
as  16. 1  percent  for  a  5-hour  forecast  using  the  hourly  updated  images.  A 
reduction  in  whole-image  RMSE  does  not  always  reflect  an  improvement  in  the 
subjective  visual  quality  of  the  forecast  images.  Indeed,  in  these  cases,  elimination 
of  bilinear  interpolation  did  not  produce  dramatically  sharper  images. 

Nevertheless,  the  decision  was  made  on  the  basis  of  objective  and  subjective 
testing  to  retain  bilinear  interpolation  in  the  cross  correlation  code. 
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The  question  of  how  many  times  to  apply  bilinear  interpolation  in  producing  a 
forecast  image  remained.  The  cross  correlation  technique  was  originally 
configured  iteratively,  such  that  a  forecast  for  hour  H  was  based  on  the  previous 
hour's  (H  minus  1)  forecast  image,  which  had  in  turn  been  based  on  the  forecast  at 
H  minus  2,  and  so  on,  resulting  in  forecast  grayshade  values  that  had  been 
bilinearly  interpolated  several  times.  Alternatively,  in  an  effort  to  minimize 
smoothing,  we  have  attempted  to  trace  all  grayshade  trajectories  (displacement 
vectors)  back  to  their  origins  in  the  most  recent  observed  image,  regardless  of  the 
forecast  interval  involved.  This  has  produced  uniformly  sharper-looking  forecast 
images  without  sacrificing  accuracy  (as  measured  by  RMSE).  We  therefore 
recommend  the  non-iterative  approach  to  the  smoothing  problem. 

2.4  Development  and  Dissipation  of  Cloud 

The  cross  correlation  technique  is  best  suited  to  producing  forecast  satellite 
imagery  in  the  cool  season  (October  -  March  in  the  Northern  Hemisphere)  over 
midlatitude  areas.  This  is  because  the  original  technique  was  designed  only  to 
advect  existing  cloud,  and  advective  processes  often  dominate  in  extratropical 
regions  from  late  autumn  through  early  spring.  It  does  not  currently  have  the 
capability  to  extrapolate  intensity  of  grayshade  associated  with  the  cloud  growth  or 
decay  implied  in  the  differences  between  the  two  initial  satellite  images. 

Apart  from  simple  advection,  cloud  grayshade  change  could  be  the  result  of 
thickening/thinning,  vertical  ascent,  or  horizontal  spreading  of  the  cloud  as  it  either 
moves  horizontally  or  remains  stationary.  Convection  is  a  major  problem  in  this 
regard  that  is  most  prevalent  during  the  warm  season  in  the  midlatitudes  and 
throughout  the  year  in  the  tropics.  Rapid  development  and  dissipation  of  cloud 
often  (though  not  always)  associated  with  convective  activity  simply  cannot  be 
resolved  by  the  cross  correlation  technique. 

Efforts  to  incorporate  the  effects  of  convection  in  the  scheme  are 
understandable  considering  its  effect  on  the  distribution  of  cloud.  However, 
attempting  to  parameterize  developmental  and  dissipative  processes  in  the 
procedure  while  simultaneously  preserving  its  ability  to  provide  forecast  images 
fairly  quickly  (and  with  a  minimum  of  outside  data  sources)  is  unrealistic. 

The  simplest  approach  is  to  calculate  the  difference  in  grayshade  for  each  pixel 
from  the  first  to  the  second  image,  and  then  to  add  that  difference  for  each  of  the 
five  forecast  images.  While  uncomplicated,  this  method  obviously  fails  when  the 
trend  established  for  the  first  two  images  is  reversed  over  the  next  several  hours. 
Even  when  this  does  not  occur,  the  increment  (or  decrement)  of  grayshade 
determined  from  the  two  observed  images  can  be  substantial  for  many  pixels.  For 
example,  the  range  of  grayshades  is  from  zero  (black)  to  256  (white).  If  a 
particular  pixel  of  interest  has  a  grayshade  of  160  units  in  the  first  image  and  180 
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in  the  second,  the  increment  is  20  units.  By  the  5th  forecast  image,  100  units 
would  be  added  to  the  grayshade  for  a  total  of 260  units.  This  would 
automatically  be  reset  to  the  displayable  maximum  of 256,  but  the  cumulative 
effect  of  many  such  pixels  would  be  splotchy  and  artificial-looking  "cloud"  in  the 
forecast  images  (the  normal  range  of  grayshades  in  our  database  images  is  roughly 
60-170).  Analogous  results  could  occur  at  the  other  end  of  the  scale  as  well.  In 
fact,  this  extreme  stretching  of  grayshades  was  all  too  common  when  applying  this 
"grayshade  trend"  approach.  There  may  be  some  promise  in  somehow  scaling  the 
increments/decrements  so  as  to  remain  in  the  normal  range  of  grayshades  for 
satellite  images,  and  we  hope  to  look  into  this  in  the  near  future. 

Our  current  work  centers  on  inferring  (rather  than  modeling)  development  and 
dissipation  of  cloud  using  cloud  climatologies,  a  strategy  most  likely  to  be 
successful  in  areas  where  cloud  cover  has  a  strong  diurnal  signal  (such  as  the 
tropics).  Two  potential  methods  for  correcting  the  advected  grayshade  value  for 
diumally-associated  growth  or  decay  were  considered  in  this  work  unit.  The  first 
can  be  thought  of  as  a  climatological  blending  approach.  It  involves  modif^ng 
standard  cross  correlation  forecasts  of  pixel  grayshade  intensities  on  the  basis  of 
hourly  cloud  climatologies.  The  other  uses  harmonic  analysis  to  discern  a  signal  in 
the  temporal  variation  of  grayshade  over  time  at  a  given  location  or  region.  The 
diagnosed  signal  can  then  be  used  to  nudge  the  standard  cross  correlation  forecasts 
of  grayshade  intensities.  Both  methods,  described  in  detail  below,  depend  on  the 
availability  of  an  hourly  cloud  analysis  archive.  Hourly  fi-actional  cloud  cover,  as 
determined  by  the  SERCAA  cloud  analysis  algorithms,  is  archived  on  a  (25  km 
resolution)  grid  basis. 

2.4. 1  CLIMATOLOGICAL  BLENDING 

The  first  step  in  using  this  method  is  to  compute  the  hourly-averaged 
grayshade  value  over  a  sample  period  (in  our  case,  ten  days)  for  each  I/16th  mesh 
location  (where  the  1/1 6th  mesh  grayshade  value  is,  in  turn,  the  average  of  the 
gayshades  of  all  pixels  in  that  l/16th  mesh  box).  Then,  denoting  the  hourly- 
averaged  grayshade  value  {x(I,J,N)},  where  I,J  are  the  indices  of  the  1/1 6th  mesh 
grid  box  that  a  given  downstream  pixel  lies  in,  and  N  is  the  hour  index  (N=l-24 
corresponds  to  UTC=00-23),  the  assigned  grayshade  value  for  the  downstream 
point  is  given  by 

-fy,.  -  (1) 


where  n  =  UTC  hour  index  of  the  forecast  time 

i,  j  =  Cartesian  indices  of  the  downstream  pbcel  location 
i',  j'  =  Cartesian  indices  of  the  upstream  location 
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This  formulation  assumes  that  an  iterative  approach  is  used  in  the  cross 
correlation  technique;  that  is,  the  grayshade  distribution  for  the  nowcast  image  is 
always  derived  from  the  previous  hour's  image.  If  the  non-iterative  approach  is 
being  used,  as  we  recommend,  n-1  in  (1)  above  should  be  replaced  with  nO, 
signifying  the  hour  index  corresponding  to  the  initial  time  (t=0,  the  second  satellite 
image  time)  of  the  nowcast. 

Application  of  climatological  blending  to  the  Eastern  Mediterranean  and 
Eastern  Asia  ROIs  produced  less  than  satisfactory  results.  The  increment  of 
grayshade  added  to  or  subtracted  from  the  forecast  greyshades  based  on  hourly 
grayshade  climatology  was  generally  far  too  small  to  simulate  cloud  growth  or 
dissipation.  While  disappointing,  this  finding  inspired  an  alternative  technique  that 
preserves  a  wider  range  of  grayshade  values  by  eliminating  the  need  to  use  hourly 
average  grayshade  values  based  on  many  images.  This  has  the  potential  to 
produce  forecast  images  with  more  dramatic  increases/decreases  in  brightness  and 
areal  coverage  characteristic  of  cloud  development  and  dissipation.  This  approach 
will  be  tested,  and  documented  in  a  follow-up  report. 

2.4.2  HARMONIC  ANALYSIS 

Given  a  suitably  long  time  series  (for  a  given  season,  since  the  temporal 
variations  are  likely  to  be  seasonally  dependent)  of  grayshade  values  for  a  given 
point,  pixel,  gridpoint,  or  region,  harmonic  analysis  can  be  used  to  discern  the 
strength  of  the  more  significant  components  of  the  temporal  variation.  If  the  first 
several  harmonics  of  the  analysis,  given  by 

Nil  9^ 

Xk(f) = W+ sin(— /0+^  cos(— /'O)  (2) 

,=i  t'  r 


where  i  =  number  of  the  harmonic  (integer  between  1  and  N/2), 

N  =  number  of  observations  in  the  time  series 
P  =  fundamental  period  of  the  sample  (total  time  period) 

[X]  =  sample  mean 

t  =  time  index  of  the  time  series  (time  increment  of  the  elements  of 
the  time  series,  in  the  same  units  as  P) 

explain  a  dominant  portion  of  the  total  variance  of  the  time  series,  then  it  may  be 
possible  to  use  a  truncated  version  of  (2)  to  estimate  the  grayshade  value  at  a  given 
location  according  to  the  historical  record  of  its  temporal  variability.  If  the 
fraction  of  the  explained  variance  is  F,  we  might  assign  a  grayshade  value  to  the 
downstream  location  given  by 

.y,  -  (I-F)X,  +  FX,  (3) 
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where 


Xa  =  the  assigned  grayshade  value  at  the  downstream  point, 

Xt  =  the  upstream  grayshade  value, 

Xh  =  the  output  of  the  truncated  harmonic  analysis  for  the 
downstream  point. 

In  this  way,  the  computed  grayshade  value  is  a  blend  of  the  trajected  and 
historic  grayshade  values,  with  the  proportions  dictated  by  the  fractional  explained 
variance  of  the  historical  record.  See  the  Appendix  for  the  computation  of  the 
above  terms. 

In  the  harmonic  analysis,  the  fiindamental  period  P  is  the  time  duration  of  the 
sample  (in  hours)  of  the  cloud  data  used  to  develop  the  method.  We  can 
determine  the  harmonic  that  corresponds  to  a  particular  period  of  oscillation  of 
interest  using  the  expression  P/i  =  pi,  where  pi  is  the  desired  period.  For  example, 
in  a  10-day  sample  of  hourly  data,  P  =  240.  The  harmonic  corresponding  to  the 
^  dumal  signal  (pi  =  24)  would  be  I  =  10.  We  can  then  compute  the  fraction  of 
variance  explained  by  the  I  =  10  harmonic  (see  Appendix)  to  determine  its 
contribution  to  the  variance  in  the  time  series.  If  we  found  that  it  was  a  major 
contribution  then  our  harmonic  truncation  I  would  have  to  be  greater  than  or  equal 
to  10. 

In  the  application  to  the  hourly  SERCAA  pixel  grayshade  values,  if  the 
"locations"  are  chosen  as  entire  SERCAA  l/16th  mesh  boxes,  the  grayshade 
values  would  be  averaged  over  the  1/1 6th  mesh  boxes  first,  and  the  time  series 
would  be  composed  of  these  area-averaged  values.  Then  the  truncated  time  series 
would  be  computed  as  described  in  the  previous  paragraph  for  each  l/16th  mesh 
box,  and  the  equation  for  that  box  is  applied  to  any  downstream  pixel  located  in 
that  box. 

Limited  testing  of  harmonic  analysis  on  images  in  one  ROI  has  produced 
inconclusive  results  to  date.  As  with  the  climate  blending  approach  described 
above,  the  technique's  dependence  on  cloud  climatologies  and  inherent  averaging 
of  grayshade  values  may  result  in  a  failure  to  produce  realistic  changes  in  the 
appearance  of  cloud  over  time  that  is  characteristic  of  cloud  growth  and 
dissipation. 

3.  DEVELOPMENT  OF  A  GRAPHIC  USER  INTERFACE  FOR  THE 
CROSS  CORRELATION  TECHNIQUE 

It  has  become  apparent  that  any  hopes  of  eventually  making  the  cross 
correlation  technique  an  operational  tool  for  weather  forecasters  (our  ultimate 
goal)  depend  just  as  heavily  on  the  development  of  a  user-friendly  graphic  user 
interface  (GUI)  as  on  generating  more  realistic  forecast  images  by  improving  the 
technique  itself  The  developmental  interface  that  currently  exists  is  adequate  for 
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research  purposes,  but  would  almost  certainly  discourage  a  forecaster  from  using 
the  tool  in  an  operational  setting.  We  have  therefore  enlisted  the  services  and 
resources  available  at  the  Computer  Science  and  the  Psychology  Departments  at 
the  University  of  Massachusetts  at  Lowell  (UML)  to  lend  image  processing, 
computer  display  and  human  factors  experts  to  the  cross  correlation  effort. 
Working  in  close  cooperation  with  PL/GPA  personnel,  the  UML  experts  will 
produce  an  operationally  enhanced  version  of  the  cross  correlation  cloud  forecast 
system.  This  will  require  three  main  components  of  work,  expanded  upon  below: 
(1)  port  the  system  to  C  and  UNIX-based  workstations  and  optimize  the  code  for 
speed;  (2)  develop  a  window-based  user  interface  in  the  X-windows/Motif 
environment;  and  (3)  create  workstation-accessible  documentation,  help  and  a 
tutorial  to  support  user  familiarization  and  training. 

3.1  Porting  and  Code  Optimization 

The  trend  in  scientific  and  engineering  computing  is  towards  UNIX-based 
workstations,  with  the  majority  of  programs  now  being  developed  in  C  (or 
sometimes  C-H-).  Thus  it  is  advisable  to  port  the  system  to  C  from  the  current 
FORTRAN  version.  This,  along  with  code  optimization,  has  in  fact  already  been 
accomplished  through  the  UML  contract  and  resulted  in  faster  run  times. 
Moreover,  advanced  workstations  will  provide  much  faster  processing,  estimated 
at  up  to  an  order  of  magnitude  greater  than  the  current  VaxStationA^S 
environment.  This  improved  performance  can  be  used  to  accomplish  greater 
speeds  using  the  current  (384  x  384)  images,  or  to  trade  off  some  of  that  speed  to 
handle  larger  images  (up  to  1280  x  1024).  It  can  also  be  assumed  that  advanced 
workstations  will  be  deployed  in  increasing  numbers  in  the  field,  thus  enhancing 
the  availability  and  utility  of  this  system  in  the  various  intended  field  settings. 

3.2  Development  of  a  User-Friendly  Interface 

The  primary  concern  is  to  make  the  system  as  user-fiiendly  as  possible.  The 
current  system  provides  a  rather  primitive  command-line  interface.  While  adequate 
for  research  and  development,  operational  forecasters  would  experience  difficulties 
working  with  it.  The  aim  is  to  take  full  advantage  of  the  multiple  window 
capabilities  of  advanced  workstations;  specifically,  to  develop  the  interface  in  an 
X-windows/Motif  environment  with  point-and-click  buttons/menus.  The  design  of 
the  interface  will  take  place  in  parallel  with  the  porting  of  the  current  system  into 
Unix/C.  Additional  features  intended  to  make  operation  of  the  system  as  easy  as 
possible  include  full  date/time  labeling  of  each  displayed  image,  forward  and 
backward  animation  capability  at  a  user-determined  speed,  display  of  latitude  and 
longitude  of  any  point  within  an  image  via  mouse  click.  When  executed  in  a  non¬ 
realtime  mode  (delayed  until  verification  images  are  available),  "on  the  fly" 
verification  of  forecast  images  will  be  available  as  well  as  graphical  display  of  areas 
of  over-and  under-forecasting  on  a  pixel  by  pixel  basis. 
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3.3  Development  of  Documentation,  Help,  and  Tutorial 

The  remaining  concern  in  operational  enhancement  is  to  provide  the  user  with 
help  in  understanding  the  basic  process  and  in  learning  how  to  use  and  interpret  the 
cross  correlation  technique.  This  requirement  will  be  met  by  creating  on-line 
written  documentation  explaining  how  to  run  the  system,  and  through  a  tutorial 
designed  to  take  the  operator  step  by  step  through  the  process. 

4.  CONCLUDING  REMARKS 

Development  and  enhancement  of  the  cross  correlation  technique  by  PL/GPA 
is  a  direct  response  to  a  documented  Air  Force  need  for  accurate  short-range 
forecasts  of  cloud  distribution.  An  added  requirement  is  the  ability  to  successfully 
use  the  technique  in  data-restricted  tactical  environments.  Initial  study  showed 
that  the  cross  correlation  technique  demonstrated  significant  promise.  Based  only 
on  two  consecutive  geostationary  satellite  images  over  a  given  location,  the 
technique  not  only  met  the  latter  requirement,  but  was  superior  to  several 
techniques  that  track  and  forecast  only  individual  cloud  features. 

While  endorsing  cross  correlation  as  a  viable  image  forecasting  technique,  the 
initial  study  also  identified  several  problems/deficiencies  that  required  further 
research.  These  included  terrain  advection,  distortion  at  image  borders,  loss  of 
visual  contrast  over  time,  and  inability  to  forecast  cloud  development  and 
dissipation.  All  but  the  last  issue  have  been  successfully  addressed  in  the  current 
study.  Efforts  continue  to  infer  realistic  cloud  development  and  dissipation  while 
simultaneously  retaining  the  efficiency  and  simplicity  of  the  current  version  of  the 
technique. 

Finally,  with  an  eye  toward  potential  operational  implementation  of  the 
technique,  work  is  currently  being  done  under  contract  to  port  the  entire  cross 
correlation  code  to  the  UNIX  platform  to  maximize  portability  and  to  produce  a 
user-fiiendly  GUI.  We  expect  that  this  effort  will  facilitate  inclusion  of  the 
technique  as  one  more  tool  for  forecasters  to  use  in  predicting  the  weather. 
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Table  1 .  Comparison  of  Mean  RMSE  and  Standard  Deviations  for  Eight 
Cross  Correlation  Cases  With  and  Without  Bilinear  Interpolation,  Using 
Observed  Satellite  Imagery  at  30  and  60  minute  Intervals.  RMSE  and  standard 
deviations  are  expressed  in  grayshade  units. 


30  Minute  Interval 

Witl 

lout  Bilinear  Ini 

terp. 

With  Bilinear  Interp. 

Fcst  Interval 

RMSE 

S.  D. 

Fcst  Interval 

RMSE 

S.D. 

0.5  h 

8.69 

1.75 

0.5  h 

8.30 

1.66 

l.Oh 

12.67 

2.46 

l.Oh 

11.66 

2.28 

1.5  h 

14.98 

2.86 

1.5  h 

13.62 

2.76 

2.0  h 

16.65 

3.81 

2.0  h 

15.54 

3.35 

2.5  h 

18.34 

3.87 

2.5  h 

16.09 

3.85 

60  Minute  Interval 

Witl 

tiout  Bilinear  Interp. 

With  Bilinear  Interp. 

Fcst  Interval 

RMSE 

S.D. 

Fcst  Interval 

RMSE 

S.D. 

l.Oh 

13.14 

2.68 

l.Oh 

12.08 

2.50 

2.0  h 

17.63 

3.64 

2.0  h 

15.81 

3.76 

3.0  h 

20.43 

4.59 

3.0  h 

18.19 

4.90 

4.0  h 

22.02 

5.41 

4.0  h 

19.14 

5.89 

5.0  h 

23.97 

6.27 

5.0  h 

20.11 

6.56 
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Appendix 

Harmonic  Analysis  Method 

The  time  series  of  N  data  values  can  be  expanded  into  a  Fourier  Series  in  the  form; 

—  NH  In  lit 

X{t')=X+  z  [4sin(^7)  +  £,cos(— /r)} 

z  =  1  ^  ^ 

tvhere 

X  is  the  estimated  value  of  the  variable  whose  time  series  is  being  analyzed; 
t  is  the  elapsed  time  in  the  time  series  (in  the  same  units  as  P) 

—  1  ^ 

X  —  Z  ^  ,  the  sample  mean  of  the  sample  values  x; 

^7  =  1  ^ 

P  the  fundamental  period,  (equal  to  the  total  temporal  period  of  the  sample); 

N 

i  is  the  harmonic  index  -  there  are  —  harmonics  in  a  data  sample  of  N; 
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,  2  ,  .  ,ln 

-4,  =-  X^[i^  sin{— "^)], 

„  1  ,  /2;r. 


except 


\  N  2p  N 

%/2=0.  %/2=]^^.2^[^7COs(— yty)] 

Then  the  fraction  of  explained  variance  by  each  harmonic  is  given  by 

Fi=(\ii)(a}+b})isI 

where  5"^  =(1 /AO  E  ,  the  sample’s  standard  deviation 

7  =  1 

Finally,  the  fraction  of  explained  variance  over  I  harmonics  is  given  by 

I  Nil 

Fj^  ZFy,  sothatF^/2=  SF;  =  1.0 
/■  =  1  /  =  1 

Usually,  I «  N/2,  since  most  of  the  explained  variance  should  be  in  the  first 
several  harmonics.  If  not,  it  is  an  irregular  oscillation. 
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